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1. INTRODUCTION 

It has long been appreciated that the simple bi rth-and-death Markov 
process often provides an adequate initial model for the behavior of 
service systems, populations, epidemics, and many other stochastic 
systems; see Feller [3] for an early classic account. Refinements, 
particularly in the modeling of service systems, have typically involved 
the replacement of exponential service times, also assumed independent, by 
independent random variables of "general" or arbitrary distribution, 
replacement of independent exponential inter-arri val times by independent 
"general" random variables, or both. Isolated examples in which arrival 
process parameters are allowed to change deterministical ly in time have 
also been studied. 

Randomly appearing fluctuations in system environment, associated 
with weather or other change in physical surroundings , personnel changes, 
alteration of system usage intensity, etc., are likely to be reflected 
in changes in observed failure and repair rates in a system reliability 
context. 

1.1 Examples : In order to illustrate the ideas summarized above, we 

introduce two specific birth and death models in random environments. 
Numerical results and discussion will be found in the last section. 

Repairman Model . Suppose any of m machines that are in use fail incependent- 
ly in Markovian fashion at rate { X ^ , t > 0), and, if failed, experience 
repair at rate {y t , t > 0). In turn, the rates {X^} and {y ^ } are them- 
selves finite-state Markov processes independent of the state of the 
machines, so that if n is the number of machines on repair, and j identifies 
the environmental state, then (n,j) is the overall state variable of the 



system; conditional upon j, n changes by one unit at a time in typical 
birth-death fashion with parameters depending on j : the probability that 
an operating machine goes down in (t,t+dt) is X-(m-n)dt + o(dt), while tne 
probability that a machine on repair becomes available is ^ min(R,n)dt + o(dt) , 
R being the number of repairmen. Such a setup describes groups of 
redundant equipments that all experience common environmental intensities 
simultaneously; the environmental changes are reflected in the numerical 
values of the failure and repair rates that prevail at any time point. 

The model described also may be used to represent the behavior of a net- 
work of m timesharing computer terminals that independently send messages or 
programs to a common central computing facility. Suppose the rate of message 
transmission changes with external activity, i.e. responds to an occasional 
period of unusual activity, perhaps a crisis situation. In tnis case the 
“constant" demand rate A (assumed equal for all terminals) switches almost 
instantaneously to a higher value, switching back to normal after the crisis 
elapses, but doing so repeatedly. Left alone, the central processor may well 
continue processing at the original rate, allowing congestion to simply build 
up, and eventually drop again when the demand lapses. On the other hand, 
remedial action may be taken. These phenomena suggest interesting and realis- 
tic questions concerning stochastic control, but these will not be 
considered in this paper. 

Mass Search Model . Let a group of m predators attempt to round up and capture 
a finite group of p prey. Predators search independently and with equal inten- 
sity and effectiveness, but the detection rate is allowed to depend upon exter- 
nal environmental conditions that cause relatively long-term changes in, say, 
visibility. If a prey is detected, it is followed until lost by the predator, 
after which moment it is susceptible again to search, detection, and active 
surveillance. While predators and prey move independently, as soon as a pre- 
dator begins following a particular prey the latter is removed from circula- 
tion and the remaining unattached predators continue search for the free prey. 
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Of interest is the long-run or stationary distribution of the number of 
prey under simultaneous surveillance, and also the first-passage time until 
a large fraction- perhaps all- of the prey are simultaneously under the eye 
of the predators. Here is a plausible model : the state of the system is 
(n,j), so that if the number of prey under surveillance, is n, and the 
state of the environment is j, then the probability that a free prey is 
detected is A .(p-n)(m-n)dt + o(dt), while the probability that a prey under 

J 

surveillance is lost is u-ndt + o(dt), where n < min(p,m) and A- and u* 

J J J 

change in accordance with independent finite-state Markov processes, as 
before. 

1.2 The analytic structure : We shall consider Markov processes 
{ ( Xt , Y t ) > t > on the state space {(n,j), 0 < n < N, 1 < j < K n h with a 

block-tridiagonal infinitesimal generator Q : 



A<°> A<°> 0 0 

M 0) A (i) A U) 0 



o 



o 



0 M< 2 ) A< 2 > a( 2 > 



0 



Q = 



( 1 . 1 ) 



0 



m (N-1) a (N-1) a (N-1) 

o m( n > a( n ) 
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Here A^, , . . , A^ are square matrices, respectively of order 

Kg, Kj, ..., K^. Their diagonal elements are strictly negative, all other 
elements are non-negative. The matrices 0 < n< N-l, and M^, 

1 < n < N, are rectangular, with appropriate dimensions ; their entries 
are non-negative. The rowsums of Q are equal to zero, therefore we 
have that 

-0 , 

M^e + A^e + A^e = 0, 1 < i < N-l, 

M< N >e + A< N >e - 0, 

where e denotes a column vector with unit elements. The variable Y t 
is to be interpreted as the state of the environment, and X as the 
state of the birth and death process, at time t . 

We assume, furthermore, that the Markov process Q is irreducible, and we 
denote by level n the set {(n,j), 1 < j < K n >, of states corresponding to 
the common value n for the first index. The structure (1.1) of Q permits 
the Markov process to move up or down by only one level at a time. It is 
in this respect analogous to the classical bi rth-and-death process, see 
Feller [3] or Karlin and Taylor [j>]. 

In the examples cited above, either X or Y changes each time the 
Markov process undergoes a change of state. In fact, this restriction is 
not part of the model, and we allow X and Y to change simultaneously. 

1.3 Existing literature: Infinite birth-and-death models in random 

environment have been studied for some time already; early results are 
found in Eisen and Tainiter [2], Purdue [11] and Yechiali [14]. More 
recently, Neuts [9], Chapter 6, has systematical ly examined a class of 
problems in which N = «*, and Q has a special repetitive structure : 

- a< 2 >= .... 

A^ = A< 2 >= 

M (2) = M (3) = .... 



This structure leads to matrix-geometric stationary probability vectors and 
efficient algorithmic procedures. 



5 . 



Finite models have been examined by Torrez [12], who suggests using 
numerical procedures designed to solve eigenvectors for band-matrices. 

Hajek has considered in [4], Section 5, a finite model with repetitive 
structure : 

A* 1 * - = ... - A< N - 2 ). 

a (D =a (2) . ... . A (N-1), 

m (2) =m (3) . ... =m (N-1) > 

K 1 = K 2 " - ' S-l ' K * 

Hajek determines the stationary probability vector in terms of two matrices 
R and R, of order K, which have to be iteratively computed. Finally, Keilson 
et al . [6] have considered finite models with the structure (1.1) and equal 
K n ‘$. They analyse the (Laplace transform of) first passage times distributions, 
from which they obtain equations for moments, and for the stationary distribu- 
tion. We cannot in this short space describe in detail the differences between 
our results and those in [4, 6, 12]; but we shall give some discussion at 
the end of Section 5. 

This paper presents an efficient computational approach to the analysis 
of birth-and-death models in a Markovian environment. The emphasis is upon 
obtaining numerical properties of both stationary distributions (in the next 
section), and first-passage time (in Sections 3 and 4). The computational 
algorithms are discussed, and numerical examples are given, in the last two 
sections. 

2. THE STATIONARY DISTRIBUTION 



In order to determine the stationary probability distribution, and moments of 
first passage times, it is useful to think of the Markov process as evolving in 
a certain manner. 
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For 0 < n N-l, define S p to be the restriction of the original process Q, 
observed during those intervals of time spent at level n, before the original 
process enters level n+1 for the first time, the state space of S p is 
{(n,j), 1 < j < K p } . Clearly all S p , 0 < n < N-l are transient Markov proces- 
ses The process is the restriction of the process Q to the states 

{(N,j), 1 < j < Kj } ; it is an ergodic Markov process. We denote by C p the 
infinitesimal generator of the process S n> 0 < n < N. 

In order to determine the matrices C n> we need the following result. 



Lemma 1 

Consider a Markov process on the state space { 1 ,2, . . . ,r ,r+l,r+2, . . . ,r+s} , 
with infinitesimal generator ?J : 



A 






.0 



A 

0 



( 2 . 1 ) 



where A is a square matrix of order r, A is a rectangular r by s 
matrix. The states r+1 to r+s are all absorbing. 

a. The states 1 to r are all transient if and only if the matrix A is non- 

singular. 

b. The (i,j)th entry of (-A ^ ) , for 1 <i,j ^ r; is the expected amount of 
time spent in the transient state j, starting from the transient state i, 
before absorption in any of the absorbing states. 



c . The (i ,k)th entry of (-A 1 a) , for l<i < r, r+l< k< r+s , is the 
probability that, starting from the transient state i, absorption occurs 

in the state k. 
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Proof . The first assertion is proved in Neats [9], Lemma 2.2.1. 

The matrix A’* has all nonpositive entries. The second assertion is 
a consequence of Neuts and Meier [10], Corollary 2. To determine the 
probability ^ of being eventually absorbed in state k, r+1 < k < r+s, 
starting from state i, 1 < i < r, we study the Markov chain embedded at 
instants immediately following a transition in the Markov process. 

The corresponding transition probability matrix P is given by 



'V 

P = 



i -a _1 a 
0 




A 



I 



where the matrix A is diagonal, with diagonal entries equal to those of A. 
It results from Kemeny and Snell [7], Theorem 3.3.7, p. 52, that 

P = [I - (I -A _ 1 A )] _1 (-A _ 1 A) = (-A'V; 
which completes the proof of Lenma 1. □ 



We may now determine the matrices C n . 
Lemma 2 



C = A (n) + M (n) (-C” 1 1 )A (n_1) . l<n <N. 
n n_ 1 

Proof , The equation (2.2) is obvious starting from level 0, the process S Q 
terminates as soon as the process Q enters level 1. Meanwhile, the transi- 
tions are governed by A {0] . Since the process Q is irreducible, the 
process S Q contains only transient states, and the matrix C Q is non- 
singular, by Lemma 1. Also, (-Cg) > 0- 



( 2 . 2 ) 

(2.3) 



a. 



We may now determine the generator of the process S^. Let Z^(x), 

t > 0, be defined as follows. Z^x) = i if the process S ] is in state 
(l,i) at time x. Assume that Z^x) = i, and Z^(x+ dx) = j H . 

Since j * i , a transition must have occurred in the process Q. Z^(x + dx) = j 
if and only if one of the following events have occured. 

a. The transition is from (l,i) to (l,j), this happens with probability 

a( dx + o(dx ) . 

■ *3 

b. The transition is from (l,i) to (0,k), for some k, ajid the process Q 
returns to (l,j), after spending an unspecified amount of time at level 0, 
before visiting any other state at level 1. This happens with probability 

dx . y. . + o (dx) , where y. a is the probability of moving from 

1 > K K,J K,J 

(0,k) to (l,j) before visiting any other state at level 1. 

If we set in (2.1) A = and A = A'^, it results from Lemma 1 that 




It is now clear that 




+ o (dx), 



1 < i* j < K 



and similarly, that 




+ o (dx) , 



1 < i < K l . 
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Thus, C 1 .A< 1 >»h( 1 )(-a(°))- 1 »< 0 1, 

. »(')* 

which proves (2.3) for n=l. 



Assume that (2.3) holds for n, we prove now that it holds for n+1. Let 
Z n+ i (t), t > 0, be equal to i if the process S n+ ^ is in state (n+l,i) 
at time t . 

If Z ^ (x) = i, then Z n+1 (i+di) = j*i if and only if one of the fol- 
lowing events occurs. 

a. There is a transition from (n+l,i) to (n+l,j), with probability 

A( n t^dx + o(dx ) . 

* *J 

b. There is a transition from (n+l,i) to (n,k) for some k, with probability 
M( n : ! )d T + o(di), and the process Q returns to (n+l,j), after spending an 
unspecified amount of time at levels 0,1,... ,n, before visiting any other 
states at level n+1, with probability y[ n |. 



Since S p records the visits made by Q at level n before reaching level 
n+1, we have by Lemma 1 that 



Y 



(n) 

k,j 



K 1 « (n) > k ,j- 



(2.4) 



It is now easy to prove that (2.3) holds for n+1, which completes the proof 
of the lemma . 

□ 



Wf ha ' nuw de iprmine the stationary probability vector i.e. the 
solution of the system P_ Q = 0, IP e_ = 1. We partition that vector as 
P^ = (Pg» 5 P^), where the subvectors have K n elements and 

correspond to the states of level n, 0 < n < N. 

Theorem 1 . 

The vectors P^ , 0 < n < N, are determined by the equations 

-N C N " °’ 

P,, = P^ +1 M (n+1) (-C^)» for 0 < n < N-l, 



(2.5) 

( 2 . 6 ) 



N 

(2.7) 

n=0 -n - 

Proof. From the structure (2.1) of Q, it results that the system 

_P Q = 0 may be decomposed into 

Pq A<°> + Pj M (1) = 0, (2.8) 

P n . A^ n-1 U P A^ + P . M^ n+1 ^ =0, 1 < n < N-l, (2.9) 

— n+1 — 1 

^-1 A (N_1) + ^ A (N) = o. (2.10) 



The equations (2.8) to (2.10) are matrix equivalents of the familiar 
stationary equations for bi rth-and-death processes. They may also be 

viewed as equations of probability balance. 
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From (2.9), = -Pj ^ (A^ )‘ 1 = Pj M^^-Cq 1 ). Then, by recurrence, 

using Lemma 2, the vectors P^ , 1 < n < N, satisfy the equations (2.5) and 
(2.6). Since the matrix is the infinitesimal generator of a finite 

irreducible Markov process, Equation (2.5) has a unique solution, up to a 
multiplicative constant, and that constant is determined by (2.7). 

□ 

This result suggests an algorithm to compute the stationary probability 
vector. 

Algorithm A . 

Al. Determine recursively the matrices C , 0 < n < N. 

A2. Solve the system tt^ = 0, tt^ e = 1 . 

A3. Compute recursively the vectors , n = N-l, ...» 0, using tt^ 

instead of . 

A4. Re-normal ize the vector IP so obtained. 

A complete analysis of this algorithm is deferred. We shall make some 
comments on it at the end of this section. Examples of numerical appli- 
cations to certain specific models are presented at the end of the paper. 

3. FIRST PASSAGE TIME TO HIGHER LEVELS 

We denote by T n the first passage time from level n-l to level n, 

and by t the first passage time from level n-l to level m > n : 
n,m 

T n m = inf (t > 0 : = m| Xq = n-l) , 1 < n < m < N, 



We define, for x > 0, 1 < n < N, 1 < i < 1 < j < K n , 



gj"](x) = p [T n < X, Y 



T +0 " j I X Q - n-l, Yq - i ] . 

n 
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lomu onin«j on .he state of the system at time h, a simple probabilistic 

argument shows tiat 



gj n Wh) = (1+A< n " 1 >h)g< n ].(x) + I A^’ 1 ) 

I j 1 > J J j I J 






n-2 



k=l 



M ( . n r 1} 

i ,k 



n-1 
h I g 
m=l 



(n-1), 

k,m 



g ^ n \ ( x ) 



(3.1) 



+ o(h), 



* denotes the Stieltjes convolution. Substracting 



\ ' 1 1 



sides of equation (3.1), dividing by h and letting 

(n). 



first order differential equation for q\ <(x) 



w linoie by ^ h :0 the Laplace-Stielti jes transform of gM'dx), 

G ^ n ' ( 4 ) = 17 e" g( n hx)] dx, and by G^(£) the matrix with entries 

1 ) j o uX i 

<>(U. 

We conclude from (3.1) that 

5 G (n) (0 =A (n ' 1 M n) (C)tA (n ' 1) + M ( n - 1 ) G ( n - 1 )(c)G (n )(c) . 



(n). 



(3.2) 



for 2 < n < N, 



.i^-ilarly 



G*' 1 ■ j = aWgWm + 



(3.3) 



s • • * c ■ be wri tten as 

G' :i) (0 = D 0 (C) A (0) , (3.4) 

g( n )(-) = D n j(C) a( n " n t f or 2 < n < N, (3.5) 

where = (d - A^)" 1 . ( 3 - 6 ) 

- (q . A (n) - ‘i (n) G (n) (0) _1 » l<n<N-l. (3.7) 
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Let < x. ■ j I X 0 . n-1, Y„ - 1] . 

and let G| n ’ m ^(c) denote the Laplace-Stiel tjes transform of g( n, . m ^(x) 
We readily obtain that 

G ( n »m) ( C ) = G (n,m-1)^ G^ m ^(c), 

= 5 G( k >(0, 

k=n 

for 1 < n < m < N, where we define G^ n,n "^(c) = I for all n . 



We easily prove from Lemma 2 and Equations (3.4) to (3.7) that 
D n (0) = - C; 1 , for 0 < n < N-1, 

hence G^fO) = - C'L for 1 < n < N, 

which is merely another representation of Equation (2.4) , since (G^ 
is the probability that, starting from (n-l,i), the process Q visits 
(n,j) before visiting any other state at level n. 

Let U (n ’ m ) = - — G (n,m ^(c)| 0 . 

H 5 " 



We have that 

u (n,m) = 
i >J 

and u( n * m ) = 
where u^ n ’ m ^ = 
We define 



E[t 



n ,m’ 




j | X 0 = n-1, Y 0 = i] , 



E[T n,ml X 0 = "- 1 - Y 0 = i] > 
u( n » m ) e . 

= U (n,n) = - — G (n) (s)L , and u (n) = U (n) e. 

3e ^ 



(3.8) 

(3.9) 

(3.10) 

(3.11) 
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' , jr*» ” 

I 

cted values of the first passage time to higher levels satisfy 
r ii - recurrence relations. 



1 « 

J 


= - C' 1 e. (3.12) 


l/ (nl = 


= - C^ 1 (e+M( n_1) u (n_1) ), for 2 < n < N, (3.13) 


1 > 1 1 

U 


) + G ( n ,m-1)( 0 ) u (m) > for 1 < n < m < N. (3.14) 



m . ■_ Theorem may be proved by differentiating equations (3.6) and 

expect to r and evaluating the derivatives -jjr & n U) at 5 = 0 
r< ly, argue from first probabilistic principles. First, equation 
'.12; Mows from Lemma 1 (b) and (2.2). Then let V be the elapsed 
toe unti the process enters a different level; that is, 



\l * Ulf I'L 

1. mm ? 


0 : X t t Xg > . By the strong Markov property. Lemma 1, and 

>i+i ‘0 = n ’ Y o = ^ 

n+1 ’ X V = n + 1 l X 0 = n ’ Y 0 = 1 ^ 

, i-l ; - n - 1 1 Xq - n , y q - i ] 

K n-1 

M e(i) + l [(-A (n) ) _1 M (n) ] ik EC T n l x q = n - 1 , Y Q = k] 
k- 1 


i ^ 


« ,n '”]| 1j E[T n+ ,|X 0 . n , r Q - j] 



(3.15) 
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Hence, 

[I - (-A^)" 1 M (n) (-C^ 1 )A (n_1) ] u (n+1) 

i 

= (.A< n ))-' e + (-A (n) )-' M (n) u (n) 



Multiplying both sides of (3.16) by A* n ^ results in the equation 



[A* 1 ^ + M (n) (-C^ 1 )A n ' 1 ] u (n+,) 



+ M (n) u (n) j. 



Equation (3.13) now follows from (2.3) and (3.17). 
Since 



E C T n.J X 0" , - 1 - Y 0-« 

= u (n ’ m -')(i) + G (n.m-1)( 0 ) y(m) (i) 



equation (3.14) also follows. 



(3.16) 



(3.17) 



Any number of moments of the first passage times may be similarly obtained. 
We merely state the following result for second moments, without proof. 



Id. 



i"’ m ' - E [ Vm I X 0 ' "- 1 - Y 0 ■ '!• 



V^ n) = v( n ’ n) 



n < m < N . We have 



v (1) ^(-C’V 1 ), 



» = 2 (- C'^) (I + M (n- 1 ) u ( n - 1 ))u( n ) 

+ C n"-\^ M (n " 1 M n " 1 ^, for 2<n<N, 



( 3 . 13 ) 

( 3 . 13 ) 



,(n,m) = v (n,m-l) + 2 u (n,m-l) u (m) + G (n,m-1) (0) y (m) 



( 3 . 20 ) 



x ^ i 



m < N. 



■ .r. s , theorems 1 and 2 show how the matrices C n> 0 < n < N, determined 
anna 2, play a central role in the determination of both the stationary 
hty distribution, and the moments of first passage times to higher 



tll> I v 
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Example: Consider the repairman model of Section (1.1). Assume there 

are just two environment states, denote by j = 1, 2 . Let the transition 
rate from environment state 2 to 1 be a , and from environment state 
1 to 2 be 3 . Then G^ n ^(c) satisfy the following system of equations. 



Mm) 



p ( 1 ) (r) = 1 n ( 1 \ .[■ 01 r(D /r\ 

G l,j U) " d 1 (l) + C £ j (1) d ] (1 ) + 5 6 2,j U) 



r ( 1 ) (r\ = \ o (2) + M gO ) (r) 

G 2,j d^lTu V J d 2 (l) + C G l,j Uj 



(3.21) 



/ n+ i \ A (m-n) 

g}j 1) (0 - , u\ x c MB 



d, (n) + e T ' "d“ln) +T 



u 1 [min(R,n)] , , / n+n 

+ ■ -L . -.- V - - - [G n U) G (n+1) (C)], i 

1 >J 



■ a r( n+1 Mrl 
+ d,(n) + ? G 2,j <«> ' 



(3.22) 



,(n*l) M . X 2 ( ™-"> , ,,, , M 2 C m1 n(R,n)] (n) (n+1) 

G 2,j (5) ' d,(n) + 5 ‘j 12 ' + d,(n) + e [G (5) G U)] 2,j 



■ 3 r (n+l) (rl 

+ d 2 (n)'“c G l,j (s) 



for n <_ m - 1 
where 

d-j (1 ) = A^m + a ; 



d 2 ( 1 ) A 2 m + 3 j 

d^(n) = A^(m-n) + M-|[min(R,n)] + a ; (3.23) 

d 2 (n) = X 2 (m-n) + p 2 Cmi n ( R,n) ] + 3 ; 
f 1 i f i = j , 



^ 4- L 
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The above equations can in principal be solved recursively. Finally, 
q( 1 » n+1 ) (c) = g ( 1 )( 5 ) g (2) (c) x...x G (n+1) U) • 

Similarly, the expected first passage times satisfy the following 

recursive equations. 

E[T il y 0 = n = dTHT* d^TT E[T i iv o " 2] 

tlT i l Y o ' 2] = d^ry + E[T d Y o = 1] 

« T n + ii Y o = '] = 3^y + d^y E [ T n + il Y o = 2 5 

y n [min(R,n)] f 1 

+ -^TO { E tT n |V 0 = 13 - ECT n+1 |V 0 = 1]} 

E t T n + l' Y 0 = = d^T77 + d^frrr ECT n + ll y 0 = 111 

y 9 [min(R,n)] r 'I 

{ E t T ni Y 0 ’ ^ * E [W Y 0 ■ 2 ^} 



(3.24) 
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4. FIRST PASSAGE TIME TO LOWER LEVELS 

There is some symmetry in the process Q, which we have not exploited 
yet. This we proceed to do now. Instead of the processes S n , 0 < n < N, 
defined earlier, we now consider the processes S n , 0 < n < N. 

For 1 < n < N, S p is the restriction of the process Q, observed during 
those intervals of time spent at level n, before the process Q moves down 
to level n-1 for the first time. All S n , 1 < n < N, are transient Markov 
processes. The process Sg is the restriction of Q observed at the lowest 
level; it is an ergodic Markov process. We denote by the infinitesimal 
generator of the process S p , 0 < n < N. 

We may carry out for the processes^ exactly the same analysis as we 
did for the processes S . We indicate below the main results for two 
reasons . 

a. The first passage times to lower levels must be analysed via the 
processes § n . 

b. We shall be able to describe the precise correspondence between our 
analysis, and the analysis of Neuts for the infinite quasi -birth-and- 
death process . 

The proof of the next lemma is omitted, since it is identical to that 
of Lemma 2 and Theorem 1. 



1C . 



Lemma 3 . 

The matrices C n , 0 < n < N, are recursively determined as follows. 



C = 

a 

C n = A (n) + A^ n ^(-C n+1 -1 ) M (n+1) , for 0 < n < N-l. 

The stationary probability vector = (Pq, P^, ..., P^) is determined by 
the equations 



Po c o ~ - ’ 



P* = Pp.i A (n - 1) (-C n - 1 ), for l<n<N, 



(4.1) 

(4.2) 



N 

I P e = 1. 
n — i n — 
n=0 



The (i,j)th entry of the matrix A^ n "^(-C~^) is equal to 

[» (n - 1) (-C- 1 )] 1 >;j - *i% 1> (- C n 1 )k.j • 



K a^' 1 ) 
n A i,k 



- (-Ai n 7 n ) z — — (-c; 1 )^. j 

1,1 k=l .(n-l) ,J 

i ,i 



(4.3) 



! he factor A^ n r^/(-A- n T^ ) is the probability that, upon leaving the state 
(n-l,i), the process Q moves to the state (n,k). The (k,j)th entry of 
(-C" 1 ) is the expected time spent by the process in state (n,j), starting 
from (n,k), before hitting any state at level n-l (by Lemma 1). Therefore, 
it results from (4.3) that [A^’^-C” 1 )] . . is equal to (-A^T 1 )) times 
the expected time spent in the state (n,j), before the first return to level 
n-l, given that the process Q starts in the state (n-l,i). This is exactly 
the interpretation of Neuts' matrix R for the infinite quasi-birth -and- 



Id. 



Let t „ denote the first passage time to level n, down from level 
m,n 3 

m+1, for 0<n<m<N-l : 



Vn = inf {t>0 : *t = n 1 *0 ' n,+ 1) , 



and let 



d v . 



(m,n)_ 



E [i 



m,n 

d! n) = d(. n > n) , 



Xg = m+1 > Yg — i ] , 






for 0 < n <m < N-l. The proof of the next theorem is omitted, since it 
is identical to that of Theorem 2. 



Theorem 3 . 

The expected values of the first passage time to lower levels satisfy the 
following recurrence relations. 



^( m ) = d( m,n )= - + A^ m+1 V m+1 ^), for 0 < m < N-2, 

d (m ’">= d( m ’ n+1) + G( m > n+1 )(g)^( n ) , for 0 < n < m < N-l, (4.6) 

where G (m,n) (0)= n (-C k+1 ) M (k+1 ^, for 0 < n < m < N-l, ( 4 . 7 ) 

k=m 



(4.4) 

(4.5) 



Observe that in the right-hand side of (4.7) » the left most matrix 
in the product is (-C^j), not (- C~jj) . 
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5. REMARKS ABOUT THE COMPUTATIONAL ALGORITHMS 

The algorithm A described in Section 2 is numerically stable 
under a large range of values for the entries of the generator Q 
of (1.1). The matrices (-C n *) have only non-negative entries. Moreover, 

since they measure the time spent at level n only, before moving to 
level n+1, they usually are of the same order of magnitude for each n, 

even for large values of N, if the entries of the matrices are of 

the same order of magnitude for each n, and similarly for the matrices 
A potential source of trouble exists when the A's are either very small or 

very large, compared to the other elements of Q. In that case, the expected 

times spent at one level before moving up are respectively very large or 
very small, and there is a risk of encountering overflow or underflow 
problems, when determining the matrices (-C~*). 

The steps A3 and A4 of Algorithm A are more delicate. We start step A3 
with a vector normalized by e = 1 . If N is very large, it is likely 

that the vector P M will be much smaller than tt m , and there is a real risk 
of running into overflow problems while performing step A3. In order to 
overcome this difficulty, we have merged the two steps A3 and A4, and 
re-normal ized the vectors each time a new subvector is determined (see 
Algorithm B in Appendix A). 

f n m) 

It results from Theorem 2 that one may compute the vectors u_' n ’ at 
the same time as one is preparing the evaluation of the stationary probabi- 
lity distributions, i.e. during step A1 of Algorithm A. In the numerical 
examples presented in the next section, we have determined the first passage 
times from level 0 to level m, for 1 <m < N. The corresponding algorithm 
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Finally, we compare the numerical efficiency of three approaches. 

As observed by Torrez [12], the matrix Q is a band-matrix, and there exist 
efficient numerical procedures to solve eigenvectors for band-matrices. 



The complexity of such procedures, for the matrix Q, is 0 (N,k ), where 
% 

K = max { K n , 0 < n < N] (Wilkinson and Reinsh [ 13] , p. 70) . The crucial 

step in our algorithm resides in the inversion of the matrices C n> 0 < n < N-l, 

and the solution of the system (2.5). 

N 3 

The corresponding complexity is 0( I K ). In Keilson et al. [6], 2 N matrices 

n=0 n 

have to be inverted, and N systems have to be solved. Again, the complexity 
N 3 

is 0( I K ) . Clearly, the three approaches have globally similar numerical 
n=0 n 

efficiency; in order to distinguish among them, one would have to determine the 
coefficients implicit in the 0(.) notation. However, it must be observed that 
the last two have the additional advantage of offering clear probabilistic 
interpretations of the computed quantities. 



6. NUMERICAL RESULTS 

6.1 Machine repairman model . We consider a system with N=5 machines and 
R=1 repairman. The system can be in K=2 environments. The system remains 
in the environment state j (j = 1,2) for an exponentially distributed random 
interval of time, with parameter a-. The failure rate of each machine is 

J 

equal to Aj in the j th environment, with Aj = 0.12, and X£ = 0.06. The repair 
rate of the repairman is y = 1 in both environments. 

Our objective is to measure how the rate of changes in the environment 
influences the system behaviour. In order to do so, we set = e/2, = B, 

and chose different values for $. Then the stationary probabilities of being 
in each environment remain constant and are given by y^ = 2/3 and ~ 1/3. 

We expect that if e is large, then the environment changes rapidly, and the 
system is only influenced by the average failure rate Xq = y,X,+ = ^*1* 
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On the other hand, if g is small, then the environment stays for long 
periods of time in the same state, and this should affect the dynamic 
behaviour of the system. 

We denote by £.|(g), 0 < i < 5, the marginal distribution of the 
number of machines on repair, for a given value of 6. 

We furthermore denote by n^(X), 0 < i < 5, the probability distribution 
for the classical machine repairman system, with constant failure rate X. 

In table I, we give the cumulative probabilities , corresponding to 

-5 -2 -1 

n(x), for x = x^, Xj and x 2 , and to £(g), for g = 10 , 10 , 10 , 1, 10, 

10 , 10 . We observe that the distributions of all the systems in a random 
environment are close to the distribution for the system with a unique, 

average rate Xq, with increasing differences when the environment changes 

-5 

more slowly. We also observe that £(10 ) = y ^h( X ^ ) + Y 2 _n(x 2 ), up t0 ^ ive 

decimal places. 

We conclude therefore that for the present model, the random environment 
has little effect on the marginal stationary distribution of the number of 
machines on repair. 

In order to measure the influence of the random environment on the 
dynamic behaviour of the system, we have computed the average time needed to 
reach the states [n machines on repair], 1 < n < 5, starting from the state 
[0 machine on repair]. We present on Figure 1 the results for n=5, which is 
interpretable as the "time to complete failure". The functions are as follows 

- 1 and 2 : fj(g) and f 2 (g), where fj(g) = Ettime to reach the state [all 
machines on repair] starting from the state [0 machine on repair and 
envi ronment j] }. 

- 3, 4 and 5 : g^), g(x 2 ) and g(x^), where g(x) = E{time to reach the state 
[all machines on repair], starting from [0 machine on repair]},for a system 
with constant failure rate i. 
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- 6 : y 1 ‘P 1 (3) + Y2^2^^’ et * ual t0 the stationary expected time to complete 
failure, starting from [0 machine on repair]. 

We clearly observe several typical ranges of values for £. If 6 > 10 , 
the environment changes so rapidly that the expected time to complete failure 
is equal to the time for a system with unique failure rate equal to Xg. 

For 10 _1 < £<10, the expected time to complete failure does not depend on the 

r- 

initial environment state. For £ < 10 °, the environment changes so slowly 
that complete failure is reached before a change of environment occurs. 

In fact, it seems that the system is almost completely decomposed in two 
different systems, one corresponding to each environment, with very slow 
migrations from one to the other. Also, we note that for £ < 10 , the sta- 

tionary expected time to complete failure looses any practical significance. 
Finally, the average failure rate Xg yields an overestimation of the time to 

_3 

complete failure for 10 <B<1, that is, for values of £ which are neither much 
lower, nor much higher than the failure rates. 

6.2 Mass search model . We define a reference model with p=15 prey and m 
predators. The system can be in K=2 environments, with parameters and 
<*2 = 2 a.. The detection and loss rates are respectively given by X^ = .001, 

X 2 = .005, = .02 and = .01. Thus, environment 2 is more favorable to the 

predators, since the detection rate is higher, and the loss rate is lower; 
however, environment 1 lasts on the average longer than environment 2. 

For this model, the rate of changes in the environment has an influence 
on the marginal stationary distribution of the number of prey under surveillance. 
In Table II, we give that distribution for five different models. Columns 1 
and 2 correspond to models with a unique environment (respectively environments 
1 and 2) . 
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olumns 3 and 4 correspond to the reference model, respectively with 

i. = 10" 4 (slow environment changes), and = 1 (fast environment changes). 

olumn 5 corresponds to a model with a unique, average environment : 

2 1 2 1 

^1 ^ — 3 — * 2 ’ = — 3 — ^*1 — 3 — w* 2 • Note that the distribution in column 

3 has two modes, respectively corresponding to the modes in columns 1 and 2. 

We have also measured how the system responds to increases in the effi- 

ncy of the predators. We have considered two ways for the predators to 

become more efficient. The first one is by becoming more numerous, the second 

is by increasing the probability of detecting a free prey, or by decreasing the 

probability of losing a prey under surveillance. We present on Figures 2 and 

3 the values of t n=5,10,15, j=l,2, where 
n ) j 



t . = E [time until n prey are under surveillance | at time 0, 

0 prey under surveillance, environment state is j], 

for the reference model with = .0001, and m equals 15 to 45. It clearly 
appears that environment 2 is more favorable for the predators. The curve for 
presents a plateau which we explain as follows. We denote by T^ j the 
time until n prey are under surveillance given that initially zero prey are 
surveillance, for a system with a unique environment, identical to 
onment j. For values of m less than 32, say, T^ ^ is so large (greater 

5 

10 on the average) that the reference model with = .0001 switches to 
environment 2 before 15 prey are under surveillance (the switch occurs on the 

4 

?rage after 10' units of time). Once the model is in environment 2, it takes 
2 3 

the average 10 to 10 units of time to have all preys under surveillance, 
occurs before the environment switches back to 1, and the total elapsed 

4 

cine is, on the average, approximately 10 . 
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On Figures 4 to 7, we present the values of t ., n=5,10,15, j = l,2, 

> J 

for models derived from the reference model. We consider = .0001, 
m equals to 15 and 20, and we multiply the rates X- by Vy" , we divide the 

J 

rates y- by VT » with y > 1, so that the probability ratio's "probability of 

J 

detecting / probability of losing" are uniformly multiplied by y in each state. 
These figures may be used in conjunction with Figures 2 and 3 to measure trade- 
off such as the following one. Suppose we start from the reference model with 
m=15, and that we double the number of predators to m=30. This will entail a 
reduction on t^ which can be measured on Figure 2. We can then measure on 
Figure 4 that the probability ratio must be multiplied by a factor y =3 in order 
to obtain the same reduction. 

In Table III, we give the approximate values of y which give the same reduction 
as doubling the number of predators, for m=15 and m=20. 
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APPENDIX A. 



Algorithm B : evaluation of the stationary probability distribution, the 

first and second moments of the time to reach level m, 1 <m < N, starting 
from level 0 . 

B1 First passage times from level 0 to higher levels . 

Bl.l Initial values. 



1 

o 

O 1 

t 

0 

1 

►— * 


(Equation (2.2)) 


u (1 > * (-C- 1 ) e; u' 1 ' 1 ' ~u<» ; 


(3.12) 


v<» * 2 (-C- 1 )i ( 1 hv (1 ' 1) » v'D ; 


(3.18) 


G (')(0) ~ (-Cj 1 ) A (0) ; G (1,1) (0) * G (1) (0); 


(3.11) 


(,<» . {-C - 1 ) 2 A<°>; - U< X > } 


(3.15) 


B1.2 Levels 2 to N. 
for n = 2 to N do 


{ (-C’^) -(A (n-1) + M (n_1) (-c^ 2 ) A (n ' 2) ) _1 ; 


(2.3) 


u<"> - (-C^Me + M* n " 1 V n ' 1) ; 


(3.13) 


H (l,nj „ ,d,n-l) + G (l,n- 1 ) (0) „(n) . 


(3.14) 


v ( n ) 4 . 2(-C^ 1 )(I + u^ 




(3.19) 



n . 



v (l’ n ) ^ v (l’ n "l) + 2 U ( 1 ’ n “ 1 ),j( n ) + Q(l’ n- 1)(Q) y( n ) ; 



G ( n )(0) (-c'^) A (n '^; 



G^ 1,n ^(0) «- G^ 1 ’ n ' 1 ^(0) G^ n ^(0); 

U (n ) <- (-C'^HI + M (n " 1) U (n " 1) )(-C‘^ 1 ) A (n-1 ); 

u (1,n) u^’H'^g^Vo)^ 1 ’ 11 " 1 ^) U (n) } 



(3.20) 

(3.11) 

(3.8) 

(3.16) 

(3.17) 



B2 Stationary distribution . 

B2.1 Initialization of the recurrence, 



{ C N - A (N) + M^K'lj) A (N_1) ; 



7r^ f- solution of i^C^|=£,Tr^e^=l} 



(2.3) 

(2.5) 



B2.2 Determination of , 0 < n < N. 



for n = N-l to 0 do 



t >- * s.., M ( " +1) (-c: 1 ); 



-n 



^fi+1 

(re-normal ize) 





N 






a «- 


I 


Ik 


e 




k=n 




k = n 


to 


N 


do 




-1 

a 


Ik } 


} 



( 2 . 6 ) 



TABLES 



Table 



Table 



Table 



I. Cumulative probabil i ties-number of machines on repair. 



II. Marginal distribution-number of prey under surveillance. 
Missing numbers are less than 5.10*^. 



III. Values of y giving the same reduction as doubling the 
number of predators. 



i 


0 


1 


2 


3 


4 


5 


n(x Q ) 


.56395 


.84593 


.95872 


.99256 


.99932 


1 


{(10 5 ) 


.56395 


.84593 


.95872 


.99256 


.99932 


1 


5(10 2 ) 


.56398 


.84589 


.95870 


.99255 


.99932 


1 


5(10) 


.56421 


.84561 


.95850 


.99249 


.99931 


1 


5(1) 


.56579 


.84397 


.95711 


.99202 


.99925 


1 


5(10 _1 ) 


.56910 


.84151 


.95424 


.99091 


.99908 


1 


5(10' 2 ) 


.57033 


.84078 


.95321 


.99047 


.99900 


1 


{do’ 5 ) 


.57050 


.84068 


.95306 


.99040 


.99899 


1 


n(Xj) 


.49516 


.79226 


.93486 


.98620 


.99852 


1 


0 (^ 2 ) 


.72118 


.93754 


.98946 


.99881 


.99993 


1 



Table I 



Cumulative probabilities-number of machines on repair 



n 


I 


II 


III 


IV 


V 


0 

1 


.00073 

.00820 




.00048 

.00543 


.00005 


.00003 


2 


.04019 




.02662 


.00061 


.00043 


3 


.11319 




.07504 


.00428 


.00342 


4 


.20375 


.00002 


.13524 


.01959 


.01724 


5 


.24653 


.00026 


.16400 


.06175 


.05835 


6 


.20544 


.00220 


.13776 


.13712 


.13601 


7 


.11886 


.01270 


.08407 


.21606 


.22012 


8 


.04755 


.05082 


.04948 


.24021 


.24628 




.01294 


.13835 


.05546 


.18523 


.18752 


10 


.00233 


.24903 


.08481 


.09610 


.09441 


11 


.00026 


.28299 


.09427 


.03198 


.03001 


12 


.00002 


.18866 


.06255 


.00633 


.00560 


13 




.06531 


.02160 


.00066 


.00054 


14 




.00933 


.00308 


.00003 


.00002 


15 




.00031 


.00010 







Table II 

Limiting distribution-number of prey under surveillance. 
Missing numbers are less than 5.10 ^ 



m 


‘5,1 


t 10,l 


t 15 , 1 


t 5,2 


t 10,2 


t 15 ,2 


15 


4.1 


4.4 


7.8 


4.9 


7.1 


> 15 


20 


4.1 


4.0 


3.9 


4.5 


5.6 


5.4 



Table III 

Values of y giving the same reduction as doubling 
the number of predators 



FIGURES 



Fig.l. Expected time to reach the state [all machines on repair], 
starting from [no machine on repair]. 

Fig. 2. Expected time to reach the states [n prey under observation], 
starting from [0 prey under observation, environment 1], 
varying number of predators. 



Fig. 3. Expected time to reach the states [n prey under observation], 
starting from [0 prey under observation, environment 2], 
varying number of predators. 



Fig. 4. Expected time to reach the states [n prey under observation], 
starting from [0 prey under observation, environment 1], 

15 predators, varying probability ratio. 



Fig. 5. Expected time to reach the states [n prey under observation], 
starting from [0 prey under observation, environment 2], 

15 predators, varying probability ratio. 



Fig. 6. Expected time to reach the states [n prey under observation], 
starting from [0 prey under observation, environment 1], 

20 predators, varying probability ratio. 



Fig. 7. 



Expected time to reach the states [n prey under observation], 
starting from [0 prey under observation, environment 2], 

20 predators, varying probability ratio. 




3-l. Expected time to reach the state [all machines on repair], startinq frt 
[no machine on repair]. 




Fig. 2. Expected time to reach the states [n prey under observation], starting 

from [0 prey under observation, environment 1], varying number of predators. 




Fig. 3. Expected time to reach the states [n prey under observation], starting 

from [0 prey under observation, environment 2], varying number of predators. 




Fig. 4. Expected time to reach the states [n prey under observation], starting 
from [0 prey under observation, environment 1], 15 predators, varying 
probability ratio. 




Fig. 5. Expected time to reach the states [n prey under observation], starting 
from [0 prey under observation, environment 2], 15 predators, varying 
probability ratio. 




Fig. 6 . Expected time to reach the states [n prey under observation], starting 
from [0 prey under observation, environment 1], 20 predators, varying 
probability ratio. 




Fig. 7. Expected time to reach the states [n prey under observation], starting 
from [0 prey under observation, environment 2], 20 predators, varying 
probability ratio. 



DISTRIBUTION LIST 



NO. OF COPIES 



Library, Code 0142 
Naval Postgraduate School 
Monterey, CA 93940 



4 



Dean of Research 
Code 01 2A 

Naval Postgraduate School 
Monterey, CA 93940 



Library, Code 55 

Naval Postgraduate School 

Monterey, CA 93940 



2 



Professor D. P. Gaver 



183 



Code 55Gv 

Naval Postgraduate School 
Monterey, CA 93940 

Defense Technical Information Center 
ATTN: DTIC-DDR 
Cameron Station 
Alexandria, Virginia 22314 



